Temperature dependence of the specific heat and the penetration depth of 
anisotropic-gap BCS superconductors for a factorizable pairing potential 
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An explicit expression for the temperature dependence of the specific heat of clean anisotropic- 
gap superconductors is derived within the weak-coupling BCS approximation. The specific heat is 
presented as a functional of the superconducting gap on the Fermi surface. The obtained formula 
interpolates between the correct low coupling jump at T c and the low temperature behavior for 
T < r c . For isotropic superconductors the formula gives a relation between the specific heat and 
the superconducting gap. For anisotropic superconductors, the interpolation formula incorporates 
averaging of powers of the gap anisotropy function over the Fermi surface and provides a suitable 
set for fitting model Hamiltonians to experimental data. The work of the interpolation formula 
is illustrated by (i) the Pokrovsky formula for the specific heat jump, (ii) Gor'kov and Melik- 
Barkhudarov formulas for the Ginzburg-Landau coefficients, (iii) the Moskalenko two-band formula 
for the specific heat jump, (iv) the temperature dependence of the specific heat for the two-band 
model, applicable to MgB2, (v) the two-dimensional d-wave model, applicable for YBa2Cu307_a, 
and (vi) the Zhitomirsky and Rice triplet p-wave model with horizontal line nodes for Sr^RuCU. 
The temperature dependence of the penetration depth is illustrated by fitting the general theoretical 
formula to the experimental data for MgB2, YBa2Cu307_^, and the triplet superconductor S^RuCU. 

PACS numbers: 74.20.De, 74.25.Bt 
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I. SPECIFIC HEAT 



Virtually all recently studied superconductors ex- 
hibit considerable anisotropy of the superconducting gap 
A p (T) over the Fermi surface e p = Ep. Despite the 
strong coupling effects and influence of disorder, which 
are all essential as a rule, for a qualitative analysis it is 
particularly useful to start with the weak-coupling BCS 
approximation for clean superconductors. In this case, 
very often model factorizable pairing potentials give an 
acceptable accuracy for the preliminary analysis of the 
experimental data. 

The aim of the present work is twofold. Firstly, we 
shall derive an explicit interpolation formula for the tem- 
perature dependence of the specific heat C(T). The for- 
mula is formally exact for factorizable pairing kernels 
which are consequence of the approximative separation in 
superconducting order parameter derived in BCS weak- 
coupling approximation by Pokrovskii. 1 Our formula re- 
produces the specific heat jump derived by Pokrovskii 1 
for arbitrary weak coupling kernels and Gor'kov and 
Melik-Barkhudarov 2 results for the Ginzburg-Landau 
(GL) coefficients of an anisotropic superconductor. That 
is why we believe that the suggested formula can be 
useful for the analysis of experimental data when only 
gap anisotropy and band structure are known. Sec- 
ondly, within the same system of notions and notation 
we present the recent results by Kogan 3 for the penetra- 
tion depth A(T), and propose for the zero-scattering case 
new formulas which may be used for experimental data 
processing. 

We begin with the entropy of a Fermi system per unit 



volume divided by the Boltzmann's constant ks 



S(T) = —2 n p \nn p + (1 — n p ) ln(l — n p ), (1) 

where the factor 2 takes into account the spin degen- 
eracy and the overline denotes integration over the D- 
dimensional momentum space 



f P = 
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■/(P). 
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The Fermi filling factors of independent Fermions 
_ 1 _ E p 



exp(2z p ) + 1 : 



2T' 



(2) 



(3) 



where T is the temperature times ks, are expressed by 
spectrum of superconductor 



E p = 



|A P | 2 , 



E, 



(4) 



Here we have to emphasize that for a model factorizable 
pairing potential V p , q oc XpXq the gap function is always 
separable as a product of a temperature dependent func- 
tion which can be associated with the GL order param- 
eter Q(T) and a rigid temperature independent function 
of the momentum \p- The nontrivial results 1 is that this 
separation of the variables is asymptotically correct in the 
BCS weak-coupling limit for an arbitrary kernel which is 
generally non factorizable. In fact, a factorizable kernel 
is a fairly unnatural property which, however, can occur 
if the pairing interaction is local, intra-atomic and lo- 
cated in a single atom in the unit cell. This is the special 
case of the s-d interaction at the copper site(s) in the 
Cu02 plane; 4 The separability ansatz, though, shall be 



2 



employed here to obtain a general interpolation formula 
formally exact for factorizable kernels. We assume that 
the gap anisotropy function x P ^ s known, either as a re- 
sult of solving the general BCS equation at T c , inferred 
from experimental data processing, or merely postulated 
within some model Hamiltonian, which is often the case 
for the high-temperature and exotic superconductors. 

With the above remarks, we will derive C(T) for the 
separable gap 



A P (T) = Q(T) Xp 



(5) 



and a factorizable kernel. 5 We apply the ansatz (5) to 
the BCS gap equation 6 



A P (T) 



and use the convention that a positive sign of V p _ q cor- 
responds to attraction of charge carriers and a negative 
potential energy of interaction. Substituting here 



V, 



p.q ~ Cr\pXq 



(7) 



and introducing i] = \Q\ 2 we obtain a transcendental 
equation for the temperature dependence of the gap Q(T) 



GA(rj,T) = l, 



(8) 



A( V ,T) 



Xp tanh z p 

2 En 



where we have used the identity 1 — 2n p = tanh z p and 
the coupling constant is defined by G = 1/A(0,T C ). De- 
tails on the derivation of the trial function approximation 
Eq. (7) and the numerical solution of Eq. (8) for T c <wb 
are given in Appendix A. 

For the specific heat of the superconducting phase per 
unit volume divided by ^ we have 



C{T) = Td T 5(?7(r),T) = 2E p d T n p = C v + C A , (9) 

where d F = d/dT. Here C v is the "normal" part of the 
specific heat 



C v {T)=T{d T S) v = —g c {z p ), 



(10) 



where 



9c{z) 



6 z 1 



7T 2 COSh 2 Z ' 



J — ( 



g c (z)dz = 1, (11) 



and (dr ■ ■ )n is the temperature differentiation for con- 
stant order parameter. For zero order parameter, r/ = 
at T c and above, C„ is just the specific heat of the normal 
phase C N {T) = C v {T,r] = 0). 
Introducing 



a(r,,T) = -(OtA^ 



<d v s) T = (12) 



where 



9a{z) = 



/oo 
g a {z)dz = 1, 
-OO 



(13) 



2 cosh 2 z ' 

the other term of the specific heat 

C A = Td n S{n,T)dTV{T) (14) 

can be written as 

C A = a(r 1 ,T)[-d TV (T)}9(T c -T). (15) 

Eq. (12) is actually a Maxwell-type equation d^drF = 
drd n F, where F is the free energy: S — —(dTF) v , A = 
-(dr,F) T ; cf. Ref. 7. 

Differentiating Eq. (8) we obtain dA — and 



-d T v(T) 



(d T A) v 



(d v A) 7 



(16) 



V(T) 



where the functions a and b represent a generalization of 
the GL coefficients for arbitrary temperature and order 
parameter 



7C(3) 



b(r,,T) = -(d n A) T = ^±2-_Xp9b(zp 
ir 2 If tanh z 1 



cosh 2 z 



f 

j — ( 



9b{z)dz 



(17) 
(18) 
(19) 



and C is the Riemann zeta function. Then 

-l 2 



C A 



4tt 2 



Xl 9a(Zp) 



b 7C(3) xpg b (z P )) 



■ 6(T C - T) (20) 



and 



12 



X 2 p9a{Zp) 



C » 7 ((Vx P g b (z P ) 9c(z P ) 



6{T C -T). (21) 



The functions gi{z p ), i — a,b,c, introduced in Refs. 7 
and 8, have sharp maximum at the Fermi surface and in 
a good approximation we have 



XX^)~27V F { X n p n (y p )), 
where 

/oo 
9iWx 2 + y 2 )dx, 
-oo 

n(0) = 1, ri(oo) = 0, 



x 



2T' 
i = a, 6, c. 



(23) 



We define averaging over the Fermi surface 



Up) 



VF 



v F = v(E F ) = 8(£ p ), 



(24) 



3 



where vp is the density of electron states per unit energy, 
volume and spin at the Fermi level. In such a way we 
obtain 



and 



CAT) = ^Tvp (r c (y p )} 



Ca 12 (X 2 p r a (y p )) 2 0(T C - T) 
C v 7C(3) (X$r b (y p )) (r c (y p )) ' 



(25) 



(26) 



At T c , where the gap is small and rj(0) = 1 this formula 
gives the Pokrovskii 1 result for the reduced specific heat 
jump 



AC 



2\2 



12 (Xp) 



C N (T C ) 7C(3) ( Xp ) ' 



(27) 



For the GL coefficient Eq. (12) and Eq. (17) the approx- 
imation (22) gives 



a( V ,T) = ^ (x 2 p r a (A P /2T)> , 
Kv,T)= 7 -^( Xp r b (A p /2T)) 



(28) 



Then the specific heat takes the simple GL form for ar- 
bitrary temperatures 



C{T)=C v {n,T)+T 



a 2 far) 

b( V ,T) 



6{T C -T). (29) 



Here, for the functions on the right-hand side we have 
substituted the thermal equilibrium value of the order 
parameter -q(T) = |<2(T)| 2 , obtained from the solution of 
Eq. (8). This BCS formula (29) is an example how good 
the physical intuition was in the phenomenology of super- 
conductivity. According to the Gorter-Casimir 9 model 
the specific heat is a sum of a "normal" part and another 
therm, governed by the temperature dependence of the 
order parameter and having exactly the GL form. The 
Gorter-Casimir two fluid model has very simple physical 
grounds. In the self-consistent approximation, the en- 
tropy S(T, A(T)) is a function of the temperature and 
a temperature dependent order parameter A(T). The 
temperature differentiation C(T) — T(dS/dT) inevitably 
gives two terms in Eq. (9) . According to the general idea 
by Landau, 10 the order parameter is an adequate notion 
for description of second order phase transitions, regard- 
less of the concrete particle dynamics. The e-expansion 
by Wilson and Fisher is only an ingenious realization of 
the same Landau's idea when the influence of fluctuations 
is essential. 

Again, at T c the general formulas Eq. (28) give the 
Gor'kov and Melik-Barkhudarov 2 result for the GL coef- 
ficients 



a (0,T c ) = ^( X 2 >, 



K0,T c )= 7 -l^f(x$). (30) 



This result can be directly derived 7 from the variational 
free energy F(r), T) of the superconductor which close to 
T c has the GL form 

F GL ( V ,T) « a(0,T c ) (T-T c )\Q\ 2 + h(0,T c )\Q\\ (31) 

The simplest method to calculate the GL coefficients is to 
differentiate 7 the free energy after a u-v transformations 
F(r],T) = H-TS. Then 



a(0,T c ) =(d v F) T (r) = 0,T = T c ), 
b(0,T c ) =(d 2 F) T ( V = 0,T = T c ). 



(32) 



If a Van Hove singularity (VHS) is close to the 
Fermi level the formulas for GL coefficients are slightly 
modified 8 



a(0,T c )=^/ v{E F + 2T c x)g a {x)dx, (33) 

-*c J-oo 

7C(3)(x' 



b(0,T c ) 



C v {Tc) 



+ 00 



v(E F + 2T c x)g b {x)dx, 



+ 00 



v(E F + 2T c x) g c (x) dx. 



Some important references on the influence of the VHS 
on the properties of superconductors, and pioneering 
works on the two-band model are given in Ref. 8. Let 
us evaluate the upper limit which can give a VHS. Let 
us take ID density of states v(E) oc 1 / \JE — £Vhs and 
Ep = £vhs = 0; there is no doubt that this mathemat- 
ical illustration is unphysical. In this case we have for 
the reduced specific heat jump AC/Cn(T c ), Eq. (27), an 
additional factor 



g a (x 2 )dx 



/ g c (x 2 )dx / g b {x 2 )dx 
Jo Jo 



2.51, x oc VE. (34) 



Although this mathematical example is not realistic, it 
can be seen that the VHS emulates qualitatively strong 
coupling corrections to the BCS theory: an enhancement 
of AC/C N {T C ) and 2A max (0)/T c . Another simulation 
of strong coupling effects can be demonstrated by sim- 
ple model density of states, corresponding to the case of 
layered cuprates 



= 1 + fcln- 



1 



|£ - £Vhs 
For illustration, we solve the equation 
tanh( v /^TA 5 (T)/2T) 



/ 



2VC 2 + A2(T) 



ko*: = g- 



(35) 



(36) 



taking uj d = 10, G = 1/2, and k = 10. The Z = 
(2A(0)/T c )/(2tt/7) vs E vhs /T c plot is given in Fig. 1. 
It can be seen that 7% enhancement corresponds to 



4 




FIG. 1: Z = (2A(0)/T c )/(2tt/7) vs E Y us/T c computed for 
the model density of states Eq. (35). Note that 7% enhance- 
ment corresponds to -Evhs = T c and the maximum enhance- 
ment is fa 9%. 



TVhs = T c . Thus, the influence of the VHS on the spe- 
cific heat is much stronger than on the A(0)/T c ratio. 

Let us also recall the general GL formula for the spe- 
cific heat jump at T c 



AC = T, 



Q 2 (0,T c ) 
6(0, T c ) ' 



(37) 



The two-band model provides probably the simplest 
possible illustration of the derived formula for the specific 
heat; for pioneering references on the two-band model 
see Ref. 8. The model is applicable with a remarkable 
accuracy 11 to MgB2 — a material which is in the limelight 
in the physics of high-T c superconductivity over the past 
years. 

For the normal specific heat we have 



C„(T) = -tt 2 Tu f [ Cl r c ( yi ) + c 2 r c (y 2 )] , (38) 



where 



2/1 



Ai 

2T : 



2/2 



A2 

2T' 



ci + c 2 



(39) 



and c\Vp and c 2 vp are the densities of states for the 2 
bands of the superconductor. Above T c or in the case of 
strong magnetic fields B > B c2 we have 



C N {T) = -n 2 u F T. 



(40) 



As pointed out earlier, within the weak-coupling BCS 
approximation Pokrovskii 1 has proved the general sepa- 
ration of the variables Eq. (5) which for a two-band su- 
perconductor results in a weakly temperature dependent 
gap ratio 5 = A\/A 2 — Xijxi- For MgB2 determina- 
tion of the two gaps has been carried out by directional 



point-contact spectroscopy in single crystals. One can 
see that for model evaluations the temperature depen- 
dence of the gap ratio could be neglected. 
For the moments of the gap we have 



Xr> r i 



Ap 
2T 



ci6 n rj(yi) + c 2 r l {y 2 ) 
( Cl S 2 + c 2 )™/ 2 : 



= a, b, c, 
(41) 

where the normalization is irrelevant in further substitu- 
tion in the GL coefficients. Finally for the second, GL- 
order-parameter term of the specific heat below the T c 
we obtain 



Ca(T) 



8tt 2 [c 1 5 2 r a (yi) + c 2 r a (y 2 )]' 

m VFT ~ 



CiS 4 r c (yi) + c 2 r c (y 2 ) 



(42) 



For the jump of the specific heat this formula reduces to 
the Moskalenko 13 result 



AC 



12 



(ciX? 



C 2X 2 J 



C N (T C ) 7C(3) 



cixi + c zxt 



(43) 



which is, in fact, a special case of the Pokrovskii 1 formula 
Eq. (27) applied to the two-band model. For application 
of the two-band model to the specific heat of MgB 2 the 
reader is referred to Ref. 14. 

The analysis of the specific heat for MgB 2 gives per- 
haps the best corroboration of the BCS results due to 
Pokrovskii 1 and Moskalenko. 13 Solving the Eliashberg 
equation and performing first-principle calculations for 
the specific heat of MgB 2 Golubov et al. [Ref. 15, Fig. 3] 
derived 65% reduction of the specific heat jump at T c . 
On the other hand, Eqs. (27) (43), using the param- 
eters from Ref. 15, gives (x 2 ) 2 / (x 4 ) = 58% reduction 
of the AC/C N (T C ) ratio. The 7% difference between 
those two estimates is in the range of the experimen- 
tal accuracy and the Eliashberg corrections to the BCS 
result is difficult to extracted. Unfortunately, the groups 
solving the Eliashberg equation have not compared their 
results to the classical results of the BCS theory for 
anisotropic superconductors 1 in order to analyze several 
percent strong-coupling corrections to the specific heat 
jump for MgB 2 . 

In the single band case c\ — 1 and Eq. (42) gives a 
simple relation between the specific heat and the BCS 
isotropic gap 



gen 

C N (T) 



r c (y) 



12 rl{y) 



7C(3) r b (y) ' 



(44) 



where y(T) — A(T)/2T. For anisotropic superconduc- 
tors, functions of the gap have to be averaged indepen- 
dently on the Fermi surface; this is the interpretation of 
the general formulas Eq. (28) and Eq. (29). Thus, we 
have the natural generalization 



gen 

C N (T) 
where y P {T) ■ 



= (r c (y P )) + 



12 (X 2 p r a (y p )) 2 



7C(3) (x 4 P n(y P )) ' 

A p (T)/2T = Xp Q(T)/2T. 



(45) 
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FIG. 2: Superconducting-to-normal specific heat ratio 
C(T)/Cn{T) vs the reduced temperature t = T/T c accord- 
ing to Eq. (45) computed for: (i) an isotropic-gap BCS su- 
perconductor (dashed line), (ii) a two-band superconductor 
ei = C2 = 1/2 with a gap ratio parameter S = 2.41 (dash- 
doted line) and (iii) 2D d-wave superconductor \p — cos 
tan ip = Px/Py (solid line). Note that for t > 0.2 two of the 
curves would be experimentally indistinguishable. 



FIG. 3: Comparison between the superconducting-to-normal 
specific heat ratio C(T)/Cn(T); the theoretical curve is com- 
puted following Ref. 14 with c\ = 0.49, 5 = 2.9 (solid line) 
and the experimental data for MgB2 are taken from Ref. 11 
(circles). The theoretical curve is convoluted with a Gaus- 
sian kernel Eq. (49), chosen to fit best the experimental data 
(At = 0.027). The experimental data 11 are digitized from 
Ref. 14, Fig. 3. 



For illustration, we now apply this general formula to 
three typical cases and the results are shown in Fig. 2: 
(i) the isotropic-gap BCS model \ p = 1, familiar from a 
number of textbooks; 16-18 (ii) the two-dimensional (2D) 
rf-wave superconductor \p — cos2</?, tan 92 = p y /p x ', and 
(iii) a two-band superconductor c\ = C2 = 1/2, for which 
the gap ratio parameter is taken to reproduce the same 
reduced specific heat jump of the d-wave superconductor 
(5 = VS ±VS = 2.41 or 0.41). 

The latter two models are often applied to analyze the 
behavior of Cu02 or MgB2 superconductors. Note also 
the qualitative difference: for a d-wave superconductor 
we have a quadratic specific heat at T -C T c , whereas 
for a two-band superconductor we have the exponential 
behavior C(T) oc exp(— A2/2T); see also Fig. 3 below. 

Consider now the low temperature behavior of the spe- 
cific heat per unit area for a 2D d-wave superconductors. 
Close to a node the gap is proportional to the momen- 
tum component along the Fermi contour A p (0) ~ vaPi- 
The corresponding superfluid velocity va is much smaller 
than the Fermi velocity vf, which parameterizes the de- 
pendence of the normal excitations energy £ p « VFPt as 
a function of the transversal to the Fermi contour mo- 
mentum component. For the ground state quasiparticlc 
spectrum we have E p m yj v A pf + v F p\ . It is convenient 
to introduce the dimensionless variables q\ — vaPi/2T 
and q2 = VFPt/^T. In terms of the latter we have for the 
element of the area in momentum space 



~{2irh) 2 



= 4 



(2T) 2 2irqdq 2EdE 



vavf (2-kK) 2 ttH 2 vavf ' 



(46) 



where q = \fq\ 



4 



E p /2T, and for axial sym- 



metric functions we can use polar coordinates; cf. Ref. 19. 
Here we have taken into account 4 nodal points. In such 
a way Eq. (10) gives 



a(T«T c 



16 T 1 



q 3 dq 



6.89- 



T 



irh 2 vavf Jo cosh 2 q " H 2 vavf ' 

(47) 

where we used 18C(3)/tt w 6.89; cf. Ref. 19, 
Eq. (2.9). This result together with Eq. (40) gives for 
the superconducting-to-normal specific heat ratio 



1.047 



T, 



H 2 up vavf 



t, 



(48) 



where t = T/T c is the reduced temperature. The pen- 
etration depth has a similar linear low temperature be- 
havior for d-wave superconductors. 

Very often fluctuations of stoichiometry and crystal de- 
fects make the theory of homogeneous crystal inapplica- 
ble close to the critical region. Let T c (r) be a weakly 
fluctuating Gaussian field of the space vector r. Hence, 
the simplest possible empirical model is to apply a Gaus- 
sian kernel to the theoretically calculated curve. Then 
for the heat capacity we have 



C(t) 



+ OO 



thcor 



(t) exp <^ - 



(t-t'f 
2{At) 2 



df 



AtV2TT 



(49) 



The philosophy of applying the convolution technique to 
all theoretical curves with singularities was advocated in 
the book by Migdal. Such an empirically smeared curve 
with At = 0.027 describes better the experimental data 
for MgB 2 close to T c ; T c At w 1.1 K, B c2 (0) = 2.5 T and 
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B c2 (0)At = 750 G. The result is depicted at Fig. 3, where 
the smeared theoretical curve is compared with the ex- 
perimental data. 11 In order to achieve a good fit of the 
theory to the experimental data we have treated c\ and 
5 as fitting parameters (cf. Rcfs. 11, 14, 21 and 22). The 
values used ci = 0.49 and 6 — 2.9 are slightly different 
from the set of parameters used latter for computing the 
penetration depth, but are still in agreement with dif- 
ferent spectroscopic evaluations. In order to reach the 
analogous quality of the fit of C(T) for cuprates we have 
to take into account simultaneously the gap anisotropy 
and the VHS in the general expressions Eq. (12) and 
Eq. (17). 

An analogous to Eq. (49) smearing of the fluctuation 
magnetization above T c reads 



M(B,T-T C ) 
= I M t hoor(B, T-T' c ) exp 



(50) 



{T' c ~ T e ) 
2{T c Atf 



dTL 



2irT r At 



However, for big fluctuations of T c we have to take 
into account the appearance of superconducting domains. 
Such a precise investigation of fluctuations in the mag- 
netization of Nb and Sn in the past led to the discovery 
of twinning plane superconductivity. For analytical GL 
results for twinning plane superconductivity see Ref. 23. 

Here we wish to emphasize that a large body of experi- 
mental data for B C 2(T) are strongly influenced by the dis- 
order. It is imperative to cut off a region of width T c At or 
B C 2(0)At close to B C 2(T C ) if we wish to determine B C 2(T) 
by extrapolation of properties from the superconducting 
phase or fluctuation behavior of the normal phase. Var- 
ious spurious curvatures of B C 2(T) have been reported 
merely as a result of disorder of the crystals. 



II. ELECTRODYNAMIC BEHAVIOR 

An analysis of the London penetration depth tensor, 
similar to that carried out by Kogan in Ref. 3, gives 



(A 2 (T)) Q/3 = — ^2v F {r d (y p )v a Vf3), a,f3 = x,y,z, 



where 



dep 
dp ' 



-i dv P 
m p = - - 



d 2 £ r , 



dp dp 2 

are the band velocity and effective mass and 



r d (y) = (y/vr) 2 ^ 



{y/^Y + [n + 



1 



-3/2 



(51) 
(52) 

(53) 



r d (y) « 7C(3)(j//tt) 2 « 1, r d (oo) = 1. 

For comparison, the conductivity tensor of the normal 
phase in r p -approximation reads 



<7 a/3 = 2v F e 2 (T p v a V(3). 



(54) 



For the penetration depths along the principal crystal 
axes we have in the two-band model 

K\T) = A- 2 (0) r d (yi) + A" 2 (0) r d (y 2 ), (55) 

where for uniaxial crystals like MgB2 there are only 4 
constants: A Xj i(0) = X y ,i(0), A x>2 (0) = Aj, ;2 (0), A z ,i(0) 
and A z ,2(0). These can be obtained from electron band 
calculations, 24 



{\-\T)U = —2v F c b r d (^->)(v a v ) b (56) 

£ ° C 6=1,2 21 



E (v 2 (o)w^), 



6=1,2 



where the band index b labels the leaf of the Fermi sur- 
face over which the averaging of the electron velocities is 
carried out. For a discussion and details see the review 
by Kogan and Bud'ko. 3 There is a natural "Eliashbcr- 
gization" of this result (cf. Refs. 15, 22, 25-27): 



a " 2^TA 2 

^ k(Ai+.if* 



(57) 



E 



27TTA 2 



« [a 2 K) 



+ u> 



7.2 



3/2 ■ 



where u n = (2n + 1) ttT are the Matsubara frequencies, 
w«,p = Z p (uj n )uj n , A p {uj n ) = Z p (uj n ) A p (w„) and Z p (uu n ) 
is the normalization factor. Analogous expressions can 
be worked out for the specific heat. 

For a heuristic consideration of the result by Kogan 3 at 
T = see Ref. 24. At T = the Fermi surface is shifted 
as a rigid object in the momentum space under the influ- 
ence of electromagnetic field. This shift of all conduction 
electrons explains why for the penetration depth the in- 
fluence of VHS is less essential than the influence on the 
heat capacity. The increase of the kinetic energy of all 
conduction electrons is actually the increase of the Gibbs 
free energy density AG = 2e 1 c2 \ 2 j 2 - At finite temper- 
atures the number of superfluid electrons is r d (A p /2T) 
times smaller. 

The penetration depths at T = can be also expressed 
by the optical masses and the Hall constant of the normal 
metal at high magnetic field 



£qC l<-o 



1 - = 2ef -^L_ 



i 

-i _ j £. 



K 



d "P m -i I dS P 
:EF (2nhr p _t p=EF (2^hfv p 



(58) 



£v<Ei 



d 3 P 
(27Tfi)3 



e p <E F 



d 3 p 
(27Tfi)3 



7 



the last equation being a consequence of the Gauss the- 

orem Ie p <E F ^P-k = £ P =£V dS P> where dS P is the el " 
ement of the Fermi surface oriented along the outward 
normal. For an extensive discussion on galvanomagnetic 
properties of normal metals and inclusion of hole pock- 
ets with volume density n./, for 1Z^} = e(n e — nh) see 
the textbook by Lifshitz and Pitaevskii 28 or the mono- 
graph by Lifshitz, Azbel and Kaganov. 29 The Bernoulli 
effect can be easily observed in almost compensated su- 
perconductors for which n e ss and the Hall constant 
is bigger. 

In the superconducting phase the Hall constant 72-oo 
can be determined by the Bernoulli potential 



Aip = -TZ C 



1 



2e c ; 



:AW; 



(59) 



generalization for the anisotropic case can be obtained 
by the obvious replacement X 2 j 2 — > jaX^pjp. Here we 
suppose that j <§; j c (T), j c being the critical current. If 
the magnetic field B is parallel to the surface of a bulk 
superconductor this formula gives 



Aip 



-Kr 



2^o 



(60) 



All charge carriers interact with the electric potential ip, 
but only the superfluid part oc rd(A p /2T) creates kinetic 
energy. The constancy of the electrochemical potential 
in the superconductor gives the change of the electric 
potential, i.e., the Bernoulli effect. For the temperature 
dependent condensation energy AG = —B 2 (T)/2p, the 
corresponding change of the electric potential is given by 



B 2 C (T) 

n 



(61) 



For complete determination of the Hall constant IZoc, 
the penetration depth A(T) and the optical mass of con- 
duction electrons in a clean superconductor [cf. Ref. 24, 
Eq. (20)], 



eA 2 (0) 

£ C 2 7?.oo 



(62) 



we have to investigate the Bernoulli effect for thin, 
rfflim A(T), and thick, da\ m ^> A(T), superconduct- 
ing films of the same material. M cp = 2m can be called 
effective mass of the Cooper pairs; this parameter can be 
significantly increased by disorder. 

For the temperature dependence of the electrochem- 
ical potential of the normal phase we have [Ref. 29, 
Eq. (12.16)] 



eAip 



6 v(E F ) 



(63) 



Close to a VHS the influence of the energy derivative of 
the density of states can be significant and measurable. 
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FIG. 4: In-plane normal fluid density 1 - A 2 (0)/A 2 (T) vs 
reduced temperature t — T/T c computed for 3 cases: (i) 
isotropic-gap BCS superconductor (dashed line), (ii) two- 
band superconductor MgB2 with parameters ci = 0.59, S = 
7.1/2.8 (dash-doted line) and (iii) 2D d-wave superconduc- 
tor (solid line). The experimental points for YBa2Cu3C>7_,5 
(squares) are digitized from Ref. 19 and the corresponding 
theoretical 2D d-wave curve is calculated according to Eq. (68) 
with renormalization factor Z = 1.4. Some experimental 
points for MgB2 (circles) are digitized from Ref. 30, Fig. 9; 
for details see the original work. 



The entropy and specific heat related to the volume 
density of the free energy of superconducting condensa- 
tion B 2 (T)/2fi can be determined by electric capaci- 
tor measurements, applying surface temperature oscilla- 
tions. For discussions of possible experimental setups see 
Ref. 24 and references therein. 

It is a matter of technical calculations to verify the 
identity 



(yM 2 E 



71=0 



(j//tt) 



+ 00 



1 

U+ 2 



dx 



-3/2 



(64) 



> 2 cosh 2 \J x 2 + j/ 2 
which transcribes into the form 

r a {y) +r d {y) = 1. 



= 1, 



(65) 



In such a way the electrodynamic behavior of a supercon- 
ductor can be expressed in terms of the functions, de- 
fined for description of its thermodynamic behavior. Us- 
ing Eqs. (65) and (51) we obtain 



1 



(X- 2 (T)U ( r «$)w) 



(A- 2 (0)) Q/3 



2T 



(66) 



Within the framework of London electrodynamics 
p n (T) = 1 - A 2 (0)/A 2 (T) is the normal fluid density, 
and ps(T) = A 2 (0)/A 2 (T) is the superfluid one, having 



8 



1.2 
1 

0.8 
0.6 
0.4 
0.2 







/ 






* 










































/ 










\ 








































\ 












\ 

\ 








































V 












\ 










































r 

a 
Y b 










\ 


\ 




/ 


/ 


\ 


\ 
































\ 










\ 
























r 




















A 










\ 

\ 






















r , 


























t 




\ 


\ 








































1 










\ 








\ 


\ 






























1 












N 








\ 


\ 
























/ 


/ 


























s 


V 



















0.0 1.0 2.0 3.0 4.0 5.0 

y 

FIG. 5: Plot of the Ti(y) functions (i — a,b,c,d). 
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FIG. 6: Plot of the gt(z) functions (i = a,b,c). 



total charge density ps^V^oo- For a two-band super- 
conductor, Eqs. (56) and (66) give for the penetration 
depth along the principal crystal axes 



Ps(T) 



A*(0) 



A»(T) 

W a ,l + W a: 2 =1 



2^ w a , b r d { ^ ), 



(vl) b 



6=1,2 



(67) 



For a set of parameters see the review by Kogan and 
Bud'ko. 3 We take S = 7.1/2.8 according to the spectro- 
scopic data; 14,31 see also the point contact spectroscopy 
data in Ref. 32. In Fig. 4 we compare our theoretical 
calculation with the experimental data for A(T) by Car- 
rington and Manzano. 30 Here we take c\ — 0.59 which 

gives W a ,l ~ Wa,2 ~ 0.5. 

The functions r^y) for i = a, b, c, d can be easily pro- 
grammed for the purposes of experimental data process- 
ing. The graphs of r,-(y) and the corresponding gi(z) 
functions are given in Fig. 5 and Fig. 6. The temperature 
dependence of the penetration depth A(T) is also pro- 
grammed for isotropic-gap, two-band and model 2D d- 
wave superconductor. In the 2D ci-wave case the theoret- 
ical result is compared with the experimental data 19 for 
YBa2Cu307_,5, which is also depicted at the figure. The 
linear dependence of 1 — A 2 (0) / A 2 (T) at low temperatures 
for YBa 2 Cu 3 7 _5 is discussed in Ref. 19, Eq. (2.10). For 
a 2D d-wave superconductor the general formula Eq. (51) 
gives 



Pb(T) 



A 2 (0) 
A 2 (T) 



271 



r d (Z^±cos2^ (68) 
2T 2ir 



where the temperature dependence of the order parame- 
ter is described in Appendix A. We are using an oversim- 
plified model for cuprate superconductivity for which are 
neglected (i) the anisotropy of the Fermi velocity vf{p) 
along the Fermi contour; (ii) higher harmonics of the 
gap function A p along the Fermi contour and (iii) the 
influence of VHS of the density of states slightly below 



the Fermi level. For comparison between Angle Resolved 
Photocmission Spectroscopy (ARPES) data and a lattice 
model for high-T c spectrum see Ref. 4, Fig. 3. 

Let us assume now that the order parameter for 
YBa2Cu307_,5 is Z-times higher than the BCS predic- 
tion. This could be due to the influence of VHS or, which 
is more important, strong coupling effects. Inserting here 
Z = 1.4 we can see that such a renormalization well de- 
scribes the temperature dependence of the penetration 
depth in the whole temperature interval. Finally, we have 
a good working BCS-like formula. In fact, significantly 
higher A max (0)/T c than BCS prediction is in agreement 
with the ARPES data. 



III. THE CASE FOR Sr 2 Ru0 4 

Our approach is aslo applicable to the triplet super- 
conductor Sr2Ru04; for a review see Ref. 33. We adopt 
the promising gap anisotropy model by Zhitomirsky and 

with 



Rice, 34 which gives E p = \W 



|A P | 2 



. 2 Px a 2 Py a 

sin — — cos -7— 

2ft 2ft 



2 Px® . 2 Pv™ 

cos — — sin -7- 

2ft 2ft 



cos 



PzC 



2ft ' 
(69) 

where p x a/H, p y a/H, p z c/H € (0,27r). For the Fermi sur- 
face we take a simple cylinder e p ps s(Vp^. + p y ) with 
radius ppa/h ps 0.93 n. Our calculations are depicted 
in Fig. 7. In this model calculation we have taken into 
account only one band responsible for superconductiv- 
ity. Although it is not a priori clear how "good" is this 
assumption, our curve reproduces the theoretical curve 
by Zhitomirsky and Rice 34 and passes close to the ex- 
perimental points by NishiZaki et al. 35 This promising 
success encouraged us to present the theoretical predic- 
tion for the penetration depth calculated from Eq. (66). 
According to the conclusions by Zhitomirsky and Rice 34 
their model with horizontal line nodes (see also Ref. 36) 
describes the experimental data better than a model with 
vertical line nodes. For illustration, in Fig. 7 we present 
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also our calculations for a simple 2D vertical line nodes 
model with gap anisotropy function 

*«*.(!£). (70) 

Similar model was studied by NishiZaki et al.; 35 see also 
Fig. 26 in the review by Mackenzie and Maeno. 33 

From aesthetic point of view our preferences are for the 
recent model for the gap anisotropy by Deguchi et al. 37 

\A p \ 2 oc sm 2 (p x a/h) + sin 2 (p y a/h). (71) 

Such type of anisotropy can be derived in the framework 
of quasi-two-dimensional exchange models for perovskite 
superconductivity of the type of the considered for Cu02 
plane in Ref. 4. The theoretical prediction corresponding 
to Eq. (71) is also illustrated in Fig. 7 together with the 
experimental data by Deguchi et al. 37 



IV. DISCUSSION AND CONCLUSIONS 

Let us discuss now the specific heat. We have shown 
that for factorizable kernels 5 the specific heat can be rep- 
resented as a sum of a "normal" component C U (T) and 
a term dependent on the order parameter Ca(T), which 
has the same form as in the GL theory. There is one detail 
that is worth focusing on: for the s-d model for high-T c 
superconductivity 4 the kernel is indeed separable because 
the contact interaction is localized in a single atom in the 
lattice unit cell. One should only substitute the spectrum 
of the superconductor at T < T c in the known expression 
for the GL coefficients from classical work of Gor'kov and 
Melik-Barkhudarov. 2 The final expression for the specific 
heat is a generalization of the result of Pokrovskii. 1 The 
derived formulas can be easily programmed for fitting the 
experimental data of anisotropic superconductors. For 
the jump of the specific heat at the critical temperature 
AC|t c = Ca(T~) general consideration has already been 
given in Ref. 7. The derived formula is not exact, but in- 
terpolates between the correct low-temperature behavior 
and the result by Pokrovskii 1 for the specific heat jump at 
T c . That is why we believe that our interpolation formula 
Eq. (29) can be useful for preliminary analysis of the ex- 
perimental data for the specific heat in superconductors; 
for experimental data processing the accuracy could be 
comparable, e.g., with the accuracy of the Debye formula 
for the phonon heat capacity. 

We illustrated our formulas for C(T) and A(T) for 
isotropic-gap BCS model and three of the best inves- 
tigated anisotropic-gap superconductors YBa2Cu307_a, 
Sr2Ru04 and MgB2. The nature of superconductivity 
for those superconductors is completely different: high- 
Tc and low-T c , phonon- and exchange-mediated, singlet 
and triplet Cooper pairs. In all those cases the de- 
rived formulas work with an acceptable accuracy; in 
some cases we have even quantitative agreement and for 




FIG. 7: Sr2Ru04. (a) Reduced order parameter for the Zhit- 
omirsky and Rice model Eq. (69) (solid line), the 2D vertical 
line nodes model Eq. (70) (dot-dashed line), and for the 2D 
model by Deguchi et al. Eq. (71) (doted line), (b) Specific 
heat ratio C(T) /Cn(T) for the Zhitomirsky and Rice model 
(solid line), the 2D vertical line nodes model (dash-dotted 
line), and for the Deguchi et al. model (doted line). The 
experimental points (circles) from Ref. 35 are digitized from 
Ref. 34, Fig. 1. (c) Normal fluid density 1 - A 2 (0)/A 2 (T) 
corresponding to the gap anisotropy models (69)-(71). The 
experimental points (circles) from Ref. 38 are digitized from 
Ref. 39, Fig. 2. We should note that the model with vertical 
line nodes predicts spontaneous breaking of the symmetry of 
the penetration depth in the afe-plane. 
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high-T c cuprates we have shown what the BCS analy- 
sis can give. We conclude that the statistical proper- 
ties of the superconductors [thermodynamic C(T) and ki- 
netic A(T)] are determined mainly by the gap anisotropy, 
irrespective of the underlying pairing mechanism, and 
the approximative weak coupling separation of variables 1 
A p (T) = Q{T)x P is an adequate approach. It is worth 
applying the derived formulas for C(T) and A(T) for ev- 
ery new superconductor. Often after the synthesis of a 
new superconductor single crystals are not available and 
only the data for heat capacity C(T) can help the the- 
ory to distinguish between different models for the gap 
anisotropy even before detailed spectroscopic investiga- 
tion is performed. 
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APPENDIX A: ORDER PARAMETER 
EQUATION FOR ANISOTROPIC-GAP 
SUPERCONDUCTORS 

Following Rcf . 1 , let us scrutinize the derivation of and 
the solution to Eq. (8). The gap anisotropy function will 
have non-zero values only in a narrow region near the 
Fermi surface 



x P = x p 0(ud - |£p|), 



(Al) 



Later, the differential volume in the momentum space can 
be separated to Fermi surface element dS and a normal 
element dp t 



d D p = d Pt dS=— dS, 
Returning to Eq. (8) we have 



Vf(p) 



ds r . 



dp 



(A2) 



G 



(2irH) D 



^-tanh(^)<?( WD -|^|) — =1, (A3) 



where § denotes integration over the Fermi surface. With 
the account of the energy cutoff ujp, the last reads 



G 



(2Trhy 



dS 



VF P Jo ' ~ 



- 1. 

(A4) 



According to Eq. (24) we have for the density of states 
vf = Wp) =7^0 / 6 ( £ - E F) d£ d ' 



(27T^) J 

i r ds 



v F 



(A5) 



~{2ttK) d J v F ' 



Similarly, the averaging over the Fermi surface can be 
represented as a surface integral 

(/(P)) = —<f (9 /(P)- ( A6 ) 

vp J {2irh) u vf 

In these notation Eq. (A4) reads 

/ 2 T tanh(VFTAf/2T) \ 1 = 1 
\ P { + Ap / Gv f A bcs' 

where Abcs = Gvp is the dimcnsionlcss BCS coupling 
constant. 

At T = T c , where A p = and E p = \£ p \, substituting 
x = C/2T we obtain 



r M tanh 
Jo x 



.,, . tanh.T , ! 
(X P ) I dx = 



Abcs 



M=^»l. (A8) 



Now the identity 



Jo 



tanhx , / 4 7 

dx = In I —M 

7T 



gives 



^ V (Xp)ABCS/ 



Analogously, at T = we have 
d£ 



4 

\ Jo 



^ 2 + a p(o), 



Abcs 



(A9) 



(A10) 



(All) 



Then taking into account that lod^> A p (0) we have 



I 




{ ^ 2 + A ^ 

(A12) 

As we will see later, it is convenient to modify the nor- 
malization of the order parameter and gap anisotropy 
function: 

X P a n J(Xpln|Xp|)\ 

Xp = , Q = Q Xav, Xav = exp < —j- )■ . 

Xav [ {Xp} J 

(A13) 

The renormalizing multiplier \ av is chosen in order for 
the renormalized gap anisotropy function to obey the re- 
lation 



(Xp lnXp) = 0. 



(A14) 
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For the two-band model this gives 



Xav Xl X 2 ' 



(A15) 



and one can easily verify that 

cixl ln|xi| + c 2 X2 In |X2 1 = 0. (A16) 
Similarly, using 

r' 2 7T 

/ cos 2 tp In I cos ip\ dip = — ln(e/4) (A17) 
Jo 8 



we obtain for a 2D <i-wave superconductor 



Xp(^) = ^ cos 2^, 



(A18) 



X P (<P) ln|Xp(v)l^ = 0. 



Using the approximation (A12) with a renormalized 
order parameter and gap anisotropy function, from 
Eq. (All) we derive 



Q(0) = 2ujd exp 



(Xp)Abcs/ ' 



(A19) 



This equation together with (A10) gives the well-known 
BCS relation for the renormalized order parameter for 
anisotropic superconductors 1 



2Q(0) 



2tt 

7 



3.53. 



(A20) 



We assume that the density of states v{E) is almost con- 
stant in the energy interval Ep ± 2T C . 

The renormalization does not change the gap A P (T) = 
QXp — Qx~p, but hi a sense Q{T) is the "true" BCS 
gap for an anisotropic superconductor. For T = the 
BCS model gives for d-wave superconductors A p (0) = 
A max cos 2ip, where 



2A max = 2?r 2 _ 4 2g 
T c 7 Ve~ 



(A21) 



However, for cuprates we have to take into account the in- 
fluence of Van Hove singularity and strong coupling cor- 
relations. As we fitted from the temperature dependence 
of the penetration depth for YBa 2 Cu 3 C>7_,5 we have 40% 
bigger gap A max = ZA^ ax S) and 2A max /T c w 6.0. In 
such a way the thermodynamic behavior is in agreement 
with the spectroscopic data. This is a good hint in fa- 
vor of the Landau-Bogoliubov quasiparticle picture ap- 
plied to high-T c cuprates. For MgB 2 taking c\ = 0.44 
and Ai(0) = 7.1 mcV and A 2 (0) = 2.8 we obtain 
Xi w 1.17 and X 2 « 0.46. Then Q(0) = Ai(0)/xi = 
A 2 (0)/x2 ~ 6.08 mcV = 70.6 K. For the critical temper- 
ature T c = 39 K we obtain 2Q(0)/T C w 3.62 which agrees 
with the BCS ratio (A20) within 3% accuracy as found 
in Ref. 14 



For arbitrary temperatures using the identity 
x , 2 



tanh ■ 



1 



1 




(A22) 



(A23) 



P l ^^[exp^^+l 



Substituting here 1/Abcs from Eq. (All) and taking into 
account the ojd 3> | A p (0)| approximation, Eq. (A12), we 
obtain the Pokrovskii equation 



exp ■ 



(x 2 p F(2y p )) 



(xl) 



VP ~ 7 Xav t ~ T 



Ap 
T 
(A24) 



where 



«(*) 



A P (T) _ Q(T) _ Q(T) 
A p (0) Q(0) Q( ) 



(A25) 



is the reduced order parameter < q < 1 as a function of 
the reduced temperature t = T/T c . In physical variables 
Pokrovskii 1 equation reads 



ln A P (T) 



Ap{0) +(4F(A P (T)/T)) p = 0. 



(A26) 



The function F(x) associated with the right-hand side of 
Eq. (A23) is defined by an integral, for which we have 
one integral and two different summation formulas, con- 
venient for small and large arguments 40 



F(x) 



du 



DC 
OO 



= 2 



V u 2 + x 2 [exp (Vw 2 + x 2 ) + l] 
du 



(A27) 



exp^coshu) + 1 
1 



= In — + 2tt V 



JX 



1 = 1 



(2i-l)7T y/x 2 + (21 - 1) 2 7T 2 



= -2j2(-l) n K (nx), 



2K (x^ 1) w W — e 



71=1 

where for large arguments we have the approximate for- 
mula 

! 1 9 225 \ 

~ 8x + I28x 2 ~ 3972a; 3 J ' 
(A28) 

Physically, here x = A/T, w = £/T and the upper inte- 
gration bound ujd/T has been replaced by 00. For this 
function we have the approximate formulas 

i^«i )s ,i n ^ + _L C (3)x 2 , (A29) 
"fx 8tt 2 

F(x > 1) M 2Kb (a). (A30) 
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0.0 0.5 



FIG. 8: Plot of the F(x) function. The approximations to 
F(x) for x <C 1 and x 3> 1 are given by Eqs. (A29) and 
(A30), respectively. 



The Euler constant is 7 = e c w 1.781072418 and 
C(3) w 1.202056903, where £ is the Riemann zeta func- 
tion. A plot of the function F(x) is shown in Fig. 8. In 
Appendix B a simple CH — h code for numerical evalua- 
tion of this function is provided. The use of the limes () 
function is optional; it increases the accuracy, but slows 
down the computation. For fast calculations we have to 
take only several terms of the expansions Eq. (A27). The 
colon (:) in Eq. (A24) represents an iterative assignment 
in which we use the initial approximation q = 1. 

The BCS order parameter equation Eq. (A24) is not 
specific for the physics of superconductivity. Recently, 
Abrikisov has derived the same equation for the tem- 
perature dependence of the amplitude of spin density 
waves in cuprates. 

For 2D (i-wave superconductors the Pokrovskii equa- 
tion (A24) reads 



In q = — 



Z7T 

J 2cos 2 (2ip)F 



■^006(2^)2) * (A31) 



The numerical solution for the squared reduced order pa- 
rameter q 2 (t) is shown in Fig. 9. The linear dependence 
near the critical temperature t = 1 corresponds to the 
GL approximation. In Fig. 10 the squared reduced order 
parameter for MgB2 (two-band model with c\ = 0.44, 
8 = 7.1/2.8) is compared with the experimental data 
from Ref. 31. 

As a last problem, let us derive the factorizable kernel 
(7) as a result from the BCS equation (6). For w D <C E F , 
Eq. (6) reads 



\(T)=£ V q>p A p J 


= Vp (v q<p A p J 



tanh(£ p /2T) dS p 



ta,n\i{Ep/2T) 

E u 



\2nK) D v F 

d£ p \ . (A32) 
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FIG. 9: Squared reduced order parameter Q 2 (T)/Q 2 (0) vs 
reduced temperature t = T/T c . For the two-band model, the 
ci and 8 parameters are chosen so as to simulating a d-wave 
Cu0 2 : ci = 1/2, 8 = 2.41. 



At T = T c [cf. Eqs. (A8)-(A10)] this formula gives 



A 9 (T C ) 



/ 2juj d 



(Vq,pA 



'P/P- 



(A33) 



Let us also mention the dimensions of the variables. Since 

d 3 p 



the integration J 



(27rR)> 



has a dimension of 1/volumc, 



then Vq tP , being a Fourier component of potential en- 
ergy, has dimension of energy x volume. For example, the 
Coulomb potential e 2 /r has dimension of energy and its 
Fourier transformation 
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ikT d 3 r 



(A34) 



has dimension of energy x volume. The same holds 
for the contact attraction in the BCS model poten- 
tial V(r) = —G6(r) having a constant Fourier compo- 
nent — G. The density of states vp has dimension of 
(energy xvolume) -1 , A p and E p have dimension of en- 
ergy and the Fermi surface averaging brackets (. . .) rep- 
resent a dimensionless operation. 

Let the dimensionless parameter Vq denotes the maxi- 
mum eigenvalue of the problem 



( V q,pXp)p =V Xq, 



(A35) 



and Xp is the corresponding eigenvector, with normal- 
ization (x p ) — 1. The comparison of Eq. (A35) and 



Eq. (A33) gives 



2^lod 



■ exp 



vpVo 



which is identified with Eq. (A10) and we obtain 

(XqVq,pXp)q,p 



G = V 



(x p )p 



(A36) 



(A37) 
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FIG. 10: Squared reduced order parameter Q 2 (T)/Q 2 (0) vs 
the reduced temperature t = T/T c for MgB2 (solid line) with 
a = 0.44, 8 = 7.1/2.8. The experimental points for MgB 2 
(circles) are digitized from Ref. 31. 



As the maximal eigenvalue is sought, one can apply in 
this case the Krilov iterations 

X [ n+l) ex (V q , pX ™) P , ((4" +1) ) 2 ) = 1, (A38) 

starting from some solution- like trial vector Xp' '■ Then 
the gap anisotropy function Xp 1S J us t the limit of the 

Krilov iterations Xp^ ■ 

For T = 0, the gap equation (A32) gives 



(A39) 



Within the weak-coupling BCS approximation, in the in- 
tegrant 



In 



2uo, 



|A P (0)| 



\Q(0)J 



Xp\ 



(A40) 



the first term is much bigger than the second one. For 
details we refer to the original work by Pokrovsky, 1 but 
roughly speaking ln[2wrj/A p (0)] s» const ^> 1. Within 
the latter approximation for A p (0) we obtain again the 
same eigenvalue problem and this constitutes the proof 
that the momentum dependence of the gap is rigid. 
Hence we derive the separation of the variables A P (T) w 
Q(T)x P - When the term ln|x P | in Eq. (A40) is small it 
can be treated perturbatively, and according to the nor- 
malization Eq. (A14) its influence diminish. The proper- 
ties of this approximative separation of the variables can 
be simulated by a factorizable kernel 



In other words, the factorizable approximation, 
Eq. (A41) and Eq. (7), works well when the influ- 
ence of smaller eigenvalues is small. 

Generally speaking, the separability ansatz is a low-T c 
approximation; T c should be much smaller than all other 
energy parameters: energy cutoff, Debye frequency for 
phonon superconductors, exchange integrals for exchange 
mediated superconductivity, the Fermi energy and the 
bandwidths. Room temperature superconductivity is not 
yet discovered, but the good message is that we have 
still a simple approximation acceptably working for all 
superconductors. For theoretical models the accuracy of 
the separable approximation can be easily probed when 
investigating the angle between the order parameter at 
different temperatures, e.g., 



arccos 



(A;(r)A p (r c )) 

y/(\A p (T)\*)(\A p (T c )\< 



or 



arccos 



A;(r)A p (r c ) 



« 1, 



< 1. 



(A43) 



(A44) 



V 9 , P = ^V n *M*<?> ^V 0XqX p, 



(A41) 



Those angles are just zero at T c and the expressions for 
the specific heat jump and the GL coefficients is correct. 
Only for T — ► some small deviations can be observed, 
but in that case one can treat Xv as a trial function in a 
variational approach. 

The performed analysis shows that the separation of 
the variables Eq. (5) due to Pokrovsky 1 and consequent 
factorizable kernel Eq. (7) are tools to apply the weak- 
coupling BCS approximation to anisotropic-gap super- 
conductors. The factorizable kernel gives a simple solu- 
tion to the gap equation, the nontrivial detail being that 
this separability can be derived by the BCS gap equa- 
tion. The factorizable kernel has also been discussed by 
Markowitz and Kadanoff 5 and employed, e.g., by Clem 42 
to investigate the effect of gap anisotropy in pure and 
superconductors with nonmagnetic impurities. Factoriz- 
able kernels are now used in many works on exotic super- 
conductors. However, in none of them is mentioned that 
the separability of the superconducting order parameter 
is an immanent property of the BCS theory. 1 The ac- 
curacy of the separable approximation is higher if the 
other eigenvalues of the pairing kernel are much smaller 
than the maximal one. This is likely to be the situation 
for the s-d model for layered cuprates, 4 where the s-d 
pairing amplitude J s d is much bigger than the phonon 
attraction and the other interatomic exchange integrals. 
In order for us to clarify this important approach to the 
theory of superconductivity, we have given here a rather 
methodical derivation of the Pokrovsky theory. 



(n) 

where V n are the eigenvalues and ^ p are the corre- 
sponding eigenvectors of the problem 



APPENDIX B: C++ CODE FOR COMPUTING 
THE F(x) FUNCTION Eq. (A27) 



(V q , P * p n) ) P = V n V q n \ |2) = 1. (A42) /* 
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* Calculates the F(x) function. By default this function 

* uses the epsilon algorithm (the limesQ function), but if 

* you need faster calculation with smaller accuracy, you 

* can set the optimized argument to false. Uses sum for 

* x < 0.32, hankel sum for x >= 0.32, the 

* F(x) ~= sqrt (2PI/x) . exp(-x) approximation for x > 13 

* and F(x) ~= log(PI / (GAM A * x) ) for x < 0.001 
*/ 

double F (double x, bool optimized = true) { 

if (x == 0) // return maximum possible double value. 
{ 

return 1.0E308; 

} 

if(x < 0) // F(x) is an even function. 
{ 

x = -x; 

} 

// For x < 0.001 we use the approximation 
// F(x) -= log(PI / (GAMA * x)). 
if(x < 0.001) 
{ 

return log(PI / (GAMA * x)); 

// For x > 13 we use the approximation 
// F(x) "= sqrt(2PI/x) .exp(-x) . 
if(x > 13.0) 
{ 

return sqrt (2 . 0*PI/x) * exp(-x) ; 

double F_limes (double) ; 
double F_fast (double) ; 

return optimized ? F_limes(x) : F_fast(x); 

} 

/* 

* Calculates the F(x) function without resorting to 

* the epsilon algorithm. 
*/ 

double F_fast (double x) { 

double sum=0.0, oldSum, k=1.0, eps, PI2=PI*PI; 

// For 0.001 <= x < 0.32 we use the sum 
// 0.5*F(x) = 0.5 * log(PI / (GAMA * x) ) + 
// 1 - PI/sqrt(x"2 + PI~2) + 

// 1/3 - PI/sqrt(x'2 + (3PI)~2) + 

// + 

// 1/(21+1) - PI/sqrt(x"2 + ( (21+DPI) "2) + 

// + 

if(x < 0.32) 
{ 

eps = 2.0E-8 / (x * sqrt(x)); 
sum = 0.5 * log(PI / (GAMA * x) ) ; 
do 
{ 

oldSum = sum; 

sum += 1.0 / k - PI / (sqrt(x*x+k*k*PI2)) ; 
k += 2.0; 

} 

while(f abs(sum - oldSum) >= eps * fabs(sum)); 



> 



return 2.0 * sum; 



int iPade, kPade; 

// Number of the partial sums. 

int N = 11; 

// Array, storing the partial sums, 
double Sn[12] ; 

// Array, storing the limesO function result, 
double An [12] ; 
double k = 1.0; 

// The lower bound should start at 0! 
Sn[0] =0.0; 

// For 0.001 <= x < 0.32 we use the sum 
// 0.5*F(x) = 0.5 * log(PI / (GAMA * x) ) + 
// 1 - PI/sqrt(x"2 + PI"2) + 

// 1/3 - PI/sqrt(x"2 + (3PI)"2) + 

// + 

// 1/(21+1) - PI/sqrt(x"2 + ( (21+DPI) "2) + 

// + 

if(x < 0.32) 
{ 

double PI2 = PI * PI; 
for(int 1=1 ; i<=N; i++) 
{ 

Sn[i] = Sn[i-1] + 

1.0/k-PI/(sqrt(x*x + k*k*PI2)); 
k += 2.0; 

} 



> 



return log (PI / (GAMA * x)) + 

2.0*limes(Sn, An, N, err, iPade, kPade) ; 



// For 0.32 <= x <= 13 we use the sum 

// 0.5*F(x) = K0(x) - K0(2x) + K0(3x) - ... 

// + (-l)~(n+l)*K0(nx) + . . . 

for(int 1=1; i<=N; i++) 

{ 

Sn[i] = Sn[i-1] + 

bessk0(k*x) - bessk0((k + 1.0)*x); 
k += 2.0; 

} 

return 2.0*limes(Sn, An, N, err, iPade, kPade) ; 



// For 0.32 <= x <= 13 we use the sum 

// 0.5*F(x) = K0(x) - K0(2x) + K0(3x) - ... 

// + (-l)-(n+l)*K0(nx) + . . . 

eps = 2.3E-6 * pow(x, 0.8) * exp(1.5 * (x - 0.3)); 

do 

{ 

oldSum = sum; 

sum += bessk0(k * x) - bessk0((k + 1.0) * x) ; 
k += 2.0; 

} 

while(f abs(sum - oldSum) >= eps * fabs(sum)); 
return 2.0 * sum; 



/* 

* Calculates the F(x) function using the epsilon algorithm. 
*/ 

double F_limes (double x) { 
double err; 



/* 

* Finds the limit of a series in the case where only 

* the first N+l terms are known. Method: The subroutine 

* operates by applying the epsilon-algorithm to the 

* sequence of partial sums of a series supplied on input. 

* For further details, please see: 

* T. Mishonov and E. Penev 

* Int. J. Mod. Phys. B 14, 3831 (2000). 
*/ 

double limes (double* S, double* A, int N, 

doubled err, int& i_pade, int& k_pade) 

{ 

int i , k = 1 ; 

double rLimes = S [N] , A_max = 0.0; 

err = fabs(S[N] - S[N-1]); 

i_pade = N; 

k_pade = 0; 

for(i=0; i<=N; i++) 

{ 

A[i] = 0.0; 

} 

while (N - 2*k + 1 >= 0) 
{ 

for(i=0; i<=N - 2*k + 1; i++) 
{ 

A[i] = (S[i+1] != S[i]) ? 

A[i+1] + 1.0 / (S[i+1] - S[i] : A[i+1]; 

} 

if(N - 2*k < 0) 
{ 

break; 

} 

for(i=0; i<=N - 2*k; i++) 
{ 

S[i] = (A[i+1] != A[i]) ? 

S[i+1] + 1.0 / (A[i+1] - A[i]) : S[i+1]; 

if (fabs(A[i]) > A_max) 
{ 

A_max = fabs(A[i]); 
rLimes = S [i] ; 
k_pade = k; 
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i_pade = i + k_pade; 

err = 1 . 0/A_max ; k++; 

if(S[i+l] == S[i]) } 

{ 

return rLimes; return rLimes; 

} } 

} 

} 
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